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Abstract 

Using Monte Carlo techniques, we study a simple model which exhibits a competition between 
superconductivity and other types of order in two dimensions. The model is a site-diluted XY 
model, in which the XY spins are mobile, and also experience a repulsive interaction extending 
to one, two, or many shells of neighbors. Depending on the strength and range of the repulsion 
and spin concentration, the spins arrange themselves into a remarkable variety of patterns at 
low temperatures T, including phase separation, checkerboard order, and straight or labyrinthine 
patterns of stripes, which sometimes show hints of nematic or smectic order. This pattern formation 
profoundly affects the superfluid density, 7. Phase separation tends to enhance 7, checkerboard 
order suppresses it, and stripe formation increases the component of 7 parallel to the stripes 
and reduces the perpendicular one. We verify that 7(T = 0) is proportional to the effective 
conductance of a random conductance network whose conductances equal the couplings of the XY 
system. Possible connections between the model and real materials, such as single high-Tc cuprate 
layers, are briefly discussed. 
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I. INTRODUCTION 



Systems with competing interactions have been known to show complex patterns of self- 
organization Examples of such systems in two-dimensions include the following: type I 



superconducting films in their intermediate state adsorbed monolayers on surfaces 
monomolecular film at air- water interfaces ^; ferrimagnetic garnet films 5]; and doped 
antiferromagnets, on which there are a number of reports of magnetic and charge order- 
ing BBat y y yy Q. Co„esp„ndi.g to tho. expe„.enta. .epo.s. the. have 
been a number of theoretical models intended to describe them. For example, numerous 
authors have studied Ising models that include a short range ferromagnetic interaction plus 
a longer-ranged antiferromagnetic interaction 19, 20, 21, 22, 23, 24|. The ori- 

gin of those interactions vary, depending on the particular system under study. In studies of 
magnetic films, for instance, the origin of the ferromagnetic term is the exchange interaction 
between spins, while the antiferromagnetic term has its origin in dipole-dipole interactions. 
In the context of hole-doped antiferromagnets, charge stripes are believed [2^ to originate 
from a short-ranged tendency of the holes to accumulate in regions of suppressed antiferro- 
magnetism, frustrated by a long-ranged Coulomb repulsion between them. Other workers 
3| have treated statics and dynamics of stripe phases in the presence of quenched disorder. 

Besides charge stripes, other forms of electronic inhomogeneities at low temperatures have 
been observed in doped antiferromagnets. For example, granular structures consistin g of 
superconducting domains separated by non-superconducting regions have been reported [22| 
in scanning tunneling microscopy (STM) studies of underdoped Bi2Sr2CaCu208+a;; spatial 
variations of the superconducting en ergy gap and of the local density of states spectrum has 



been observed in STM experiments 



28l | on optimally doped Bi2Sr2CaCu208+a;; and studies 



of lightly doped Ca2_xNa2;Cu02Cl2 show that electronic states within certain energy ranges 
show spatial modulations in the form of checkerboards 29^ . 

This diversity of electronic structures observed in doped antiferromagnets have prompted 
us to study an XY model with annealed disorder, where different kind of geometrical orders 
might occur, and to try to determine the interplay between those orders and superconductiv- 
ity. Specifically, we have carried out a Monte Carlo study of a site-diluted, two-dimensional 
XY-model on a square lattice, in which, besides the XY coupling, there is an additional inter- 
action between mobile spins. We have considered two types of such additional interactions: 
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(i) a screened Coulomb repulsive interaction between spins, and (ii) a repulsive interaction 
between either nearest-neighbor spins or second nearest-neighbor spins. We assume that this 
Hamiltonian governs an annealed system, in the sense that the spins are free to move under 
the influence of the repulsive interaction. They are also free to order in response to the XY 
interaction. XY models have been used to study systems such as granular superconduc- 



tors 



3G 



and the superfluid transition of helium in porous media 34]. However, to the best of our 
knowledge, there are no studies of XY models with annealed disorder as the one we present 
here. 

Our main purpose in the present work is to investigate (i) what are the low temperature 
geometrical orders that occur in an XY model with annealed disorder and different types 
of repulsion between spins, and (ii) how this rearrangement of spins affects the helicity 
modulus 7 of the system. First, we calculate the XY transition temperature for various 
choices of the parameters. This transition appears to be of the Kosterlitz-Thouless variety, 
though we have not carried out detailed flnite-size scaling tests of this hypothesis. We also 
calculate the behavior of the helicity modulus 7, which behaves like a tensor for some of the 
low-temperature phases. Finally, we determine the temperatures at which other transitions 
occur. These other transitions are induced by the Ising-like spin-spin repulsion mentioned 
above, and are of two types: (i) an "antiferromagnetic" transition into a checkerboard-like 
structure for nearest-neighbor repulsion, and (ii) a transition into a stripe phase for second- 
nearest-neighbor repulsion. In both cases, we flnd that the transition has a strong effect on 
7, which is described in detail below. 

Our calculations are carried out using standard Monte Carlo (MC) simulations. In 
some cases, we calculate critical values of the parameters that determine the type of low- 
temperature structures. We also compare our results to previous studies of site diluted 
XY- models with quenched disorder [^. The behavior in the quenched and annealed cases 
are notably different, as discussed further below. 

The remainder of this paper is organized as follows. In Section II, we describe the studied 
models. Following this, we give a brief description of our computational method in Section 
III. Our results are presented in Section IV, and a discussion of them is given in Section V. 
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II. MODEL 



We consider two distinct model Hamiltonians. The first model, which we denote Model I, 
consists of a two-dimensional (2D) square lattice with a fraction p of sites occupied by planar 
spins, and the remaining 1 —p vacant. The spins have fixed length but are characterized by 
an angle 9i, and are assumed to be described by the following classical Hamiltonian: 

H = -Y: J^, cos(^. - e,) + n,n,^^^Pt^^ (1) 

<i,j> ijtj 

Here Jij = J > if the sites i and j are both occupied, and Jy = otherwise. The first sum 
is taken over distinct pairs of nearest-neighbors < i,j >, while the second sum is taken over 
all distinct pairs of lattice sites separated by a distance < Vcut, where Vcut is some 
cutoff radius, usually taken to be a few times the screening length Vc- C is a non-negative 
constant that specifies the repulsion strength, and = 1 or is the number of spins at site 
i. The first term thus corresponds to the standard XY Hamiltonian in 2D, while the second 
term represents a screened Coulomb repulsion between the spins, which is independent of 
the angles 6i and 6j. Thus, each site has two degrees of freedom: the spin variable 6i, and 
the occupation number n^. 

Physically, each site which is occupied by a spin can be interpreted as a mobile, charged 
superconducting domain, and each vacant site as a negatively charged non-superconducting 
region. The angles 6i of the "spins" represent the phases of the superconducting order 
parameters in the i*^ domain. When two superconducting domains are close to each other, 
they couple via Josephson tunneling, which is the origin of the first term in Eq. (^. 

The second term in Eq. originates in the Coulomb repulsion between charged super- 
conducting domains. While this term might appear to violate overall charge neutrality, it 
is, in fact, consistent with that requirement, as we now show. The screened Coulomb term 
[the second term in Eq. (^] should really be written as 

V- exp(-rij/rj 

J^Coul — Z^QiQj 5 

where the sums run over all pairs of lattice sites, and charge neutrality requires that J2i li = 0. 
We assume that qi = qs or qj on the superconducting or insulating sites, where pqs + (1 — 
p)qi = O.Now note that 

Hcoui = 2^i<li- qi){qj - Qi) ^ + Qi 2^{qi + Qj) ^ 
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But in this last equation, the second term vanishes because J2i Qi = 0, while the third is just 
an additive constant. The summand of the first term is non- vanishing only on the S sites. 
Thus, Hcoui can be rewritten, to within a constant, as 

Hcoul = 2^ ninj[qs - Qi) . 

This form is indeed equivalent to the second term in Eq. with C = {qs — qiY- 

At high temperatures T, the system described by Eq. ((H) is expected to be phase- 
incoherent, and the spins will point in random directions. As T is reduced, neighboring spins 
tend to align with each other in order to minimize the system energy. This short-ranged 
attraction between spins competes with the longer-ranged screened Coulomb repulsion rep- 
resented by the second term in the Hamiltonian Eq. (0). Thus, this Hamiltonian is a simple 
model of a system with competing interactions. 

We have also studied a model Hamiltonian ("Model 11") consisting of an XY term plus 
a short-range repulsion: 

H = - ^ Jij cos{6i -6j) + A riiTij + B ^ riirij. (2) 

<i,j> <hj> <i,j>' 

Here the first and second sums are carried out over distinct pairs of nearest-neighbors 
while the third sum runs over all distinct pairs of second-nearest-neighbors. A and B are 
non-negative constants specifying the strength of the nearest-neighbor and second-nearest- 
neighbor repulsion. As will be shown, Hamiltonian (j21) at low T produces patterns similar 
to those of Model I. 

We obtain the equilibrium thermodynamics of these two models by treating (^J and 
(121) classically. At a given T, the 6'j's and n^'s arrange themselves so as to minimize the 
Helmholtz free energy F, subject to the constraint that = where p is assumed 

to be specified by the experimental conditions. Since both 6i and rii are free to change, these 
Hamiltonians describe systems with annealed disorder. By contrast, in a typical system 
with quenched disorder, the n^'s are assumed fixed, and only the ^j's are free to change. 
In this quenched case, the last two terms in the Hamiltonian would have the same value 
for any configuration of the XY spins, and thus play no role in determining the system 
thermodynamics. 



III. COMPUTATIONAL METHOD 



A. Monte Carlo Approach 

We have studied Models I and II for several different values of the parameters A, B, and 
C, and for different spin concentrations p, using a Monte Carlo (MC) approach. We used a 
square lattice, generally of size 30 x 30 sites, with periodic boundary conditions, using the 
standard MC Metropolis algorithm 3^ . In all of the studied systems we have set Jij = 1 if 
both sites i and j are occupied, and Jij = otherwise. For Model I, we considered annealed 
systems with C > 0, while for Model II, we studied two classes of parameters: (i) a nearest- 
neighbor repulsion only {A > 0; B = 0), and (ii) second-nearest-neighbor repulsion only 
{A = 0; B> 0). 

We have also reproduced some previous MC studies of 2D systems with quenched disor- 
der in order to compare with our present annealed results. For the quenched systems 
(in which the spins are distributed randomly at fixed locations on the lattice), we averaged 
over 20 different quenched disorder realizations, each with a spin concentration p. For each 
disorder realization, we started the MC run with the spins arranged in a random configura- 
tion of {^i}'s and the system at T = 1.2 (in units of J/ks)- The system was then cooled in 
steps of AT = 0.05, down to T = 0.1, and AT = 0.025 for T < 0.1. For each T, we carried 
out 5 X 10"^ sweeps through the entire lattice, taking averages over the last 2 x 10"^ sweeps, 
where each sweep consisted of an MC attempt to vary the angle of each spin. 

For systems with annealed disorder, we started the system with the spins in randomly 
chosen sites on the lattice, with randomly chosen angles 6i, at a starting T = max[2C, 1.2] 
(here C is in units of J/kB) for Hamiltonian (^, and T = 1.2 for Hamiltonian (0)- The 
system was then cooled down to a T = 0.025, in steps of AT. For Model I, we took 
AT = O.IC for T > 1.0, AT = 0.05 for 0.1 < T < 1, and AT = 0.025 for T < 0.1. For 
Model II, we used the same AT's as in the quenched case. For each T, we carried out a series 
of MC sweeps over both spin angles and spin positions. During a sweep over the spin angles, 
MC changes in the spin angles were attempted, while on sweeps over spin positions, the MC 
step consisted of an attempt to move each of the vacancies to a randomly chosen occupied 
site I35I. After each sweep over spin position, ten sweeps over spin angles were carried 
out. To find the low-T spin configuration of Model I, we took total of 2 x 10® sweeps over 
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spin angles. But to compute thermal averages for Model II, we carried out 1 x 10^ sweeps 
over spin angles, discarding the first 5 x 10^ in order to allow the system to equilibrate. 
Following this relaxation, the MC thermal averages were carried out every tenth sweep over 
spin angles, immediately before a sweep over spin positions was carried out. 

For Model I, the second sum was taken only over those distinct pairs of lattice sites which 
are separated by rij < min[3rc, L/2], where L is the linear dimension of the square MC cell. 
[Thus, the cutoff radius Vcut discussed below eq. (1) is the lesser of Svc and L/2.] This 
choice of < L/2 is necessary to avoid double counting of pairs in when periodic boundary 
conditions are used, as in the present calculation. 

For both Model I and Model II, we have calculated several equilibrium quantities. To 
characterize phase coherence, we computed the diagonal elements •jaa = x, y) of the 
helicity modulus tensor 7. These diagonal elements are the spin-wave stiffness constants 
of the XY spin system, and, in a superconductor, are proportional to elements of the su- 
perfluid density tensor. They are defined as appropriate equihbrium averages over the spin 
configuration, --^xxi for example, is defined 0| by 

<i,j> ^ <i,j> 

+ .r} A - sin(^i - (3) 

where Xi is the x coordinate of the spin on the lattice site i, is the total number of lattice 
sites, and () denotes a canonical average. 7^,^ is defined by the analogous expression with 
Xi replaced by yi. In our computations, we have taken the lattice constant a = 1. To 
characterize the spins patterns, we have also calculated 

SMT)=<\nm'>. (4) 

and 

^2(g,T) =< |n(g-)|2 > -I < n{q) > |^ (5) 

where 

AT 

^io) = J2 exp(-zg ■ fj-), (6) 

where fj is the position of the jth lattice site. For q = nx/a and q = iry/a, Si{q,T) 
and S2{q,T) probe stripe formation in the y and x directions, respectively; for q = 7t{x + 
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y)/a, they probe the formation of checkerboard patterns, and for q = 2tix/ {NxO) and q = 
2TTy/{Nya) (with A^^^ and A^j, the number of lattice sites in the x and y-directions), they are 
sensitive to phase separation of the system into large domains of occupied and vacant sites. 



B. Special Approach for Low- Temperature Helicity Modulus 

At low T, 'jxxiT) and lyyiT) can also be calculated by using a mapping between these 
quantities and the effective conductances Qe^xx and Qe^yy of a related conductance network 2^ . 
Specifically, the hehcity modulus 'jaaiT 0) in the a direction satisfies the following 
relation: 

7o<,(r = 0,p) _ ge,aa{p) 



-faa{T = 0,p=l) ge,aa{P=^)' 

Here jaaiT = 0,p) denotes the a, a component of the helicity modulus tensor at T = for a 
diluted arrangement of XY spins of concentration p and the actual T = configuration, and 
1aa(T = 0,p = 1) is the (isotropic) helicity modulus of the corresponding array at p = 1. 
9e,aa{p) ^'^^ 9e,aa{p = 1) are the conductances of a certain conductance network associated 
with the original array of diluted XY spins, and is constructed as follows: We associate with 
each spin an electrical node, and with each coupling constant Jij connecting spins at sites 
i and j we associate a "conductance" gij = Jij. The effective conductance of this network 
in the a direction {a = x or y) is denoted ge,aa- Eq- allows us to calculate the ratios of 
helicity moduli at two different concentrations p from the corresponding conductances. Our 
method of calculating the (/e's needed for this mapping is explained below. 



IV. RESULTS AND DISCUSSION 



We have carried out extensive simulations for Models I and II, considering systems of 
lattice size 30x30, and a range of A, B, C, r^, Vcut, and p, as described below. We first 
consider the low-T (T = 0.025) spin configurations for systems described by Model I. Next, 
we present a more detailed study of Model II results over a range of T. While the two 
models differ in how the repulsive interaction is truncated, several features of the low-T spin 
configurations prove to be model- independent. 
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A. Model I 



We have studied Model I [Hamiltonian (^J] for several values of the screening radius Tc- 
Low-T configurations were obtained by annealing the system starting from T = max[2C, 1.2] 
to T = 0.025, as explained in the previous section. Fig. ^ shows the T = 0.025 spin 
configurations for Tc = oo. This corresponds to a unscreened Coulomb repulsion which is, 
however, truncated at r = L/2. In this and all later Figures, the white and black squares in 
the lattice correspond to occupied and vacant sites. At all p, in the absence of a Coulomb 
repulsion (C = 0), the system phase-separates into regions of occupied and vacant sites. This 
behavior is expected, since this configuration maximizes the number of nearest-neighbors 
for each spin, and hence minimizes the internal energy. 

At C = 1, phase separation no longer occurs; instead, the spins arrange themselves 
into long, unidirectional stripes, whose width increases with increasing spin concentration 
p. For C = 1 and p = 0.3 and 0.7, we see a kind of smectic pattern in the stripes: the 
stripes seem to be arranged into layers which have a characteristic thickness, though the 
spacing between stripes is not perfectly ordered. Such smectic states have been postulated 
(in the context of a quantum model) for two-dimensional metallic stripe phases in doped 



Mott insulators 



37l | . This spin arrangement is a compromise between the clustering induced 



by the short-range attractive XY interaction and the long-range repulsion produced by the 
Coulomb interaction. As C is increased further, the system undergoes a characteristic series 
of morphology changes, from long stripes, to shorter stripe-like patterns, and eventually to 
a checkerboard pattern. At some values of p and large C, there are suggestions that the 
occupied or vacant sites have a nematic order - that is, they are arranged in short stripe-like 
patterns with a preferred direction. At certain values of p and C, we also see a pattern of 
diagonal stripes (at p = 0.7, C = 4 in Fig. 1). The checkerboard patterns are the state of 
minimum energy when the Coulomb repulsion is much stronger than the XY attraction. 

We have also sampled the low-T p-C phase diagram of Model I for finite Tc (screened 
Coulomb repulsion). Fig. |21 shows snapshots of the spin configurations for Vc = 7 and p = 0.5, 
0.6 and 0.8. The destruction of the phase-separated case at C = (shown only in Fig. [T]) by 
the screened Coulomb repulsion proceeds first by formation of long, elongated domains of the 
minority component (which, for the cases shown in Fig. |21 are vacancies). These domains 
can be seen in Fig. |2l for C = 0.1. As C is increased, those elongated domains become 
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unidirectional stripes [as in {p = 0.5, C = 1.0), (p = 0.5, C = 3.0), and {p = 0.7, C = 1.0)], 
or long but tortuous stripes which coexist with small blobs [{p = 0.8, A = 1.0)]. As in Fig.^ 
a further increase in C causes the stripes to break up and leads, for the largest values of C, 
to checkerboard patterns, of which a clear example is {p = 0.5, A = 7.0) in Fig. |21 

Figs. 121 and 0] show spin configurations analogous to those of Figs. Q and |21 Once again, 
there is a characteristic sequence of changes with increasing Coulomb repulsion, from phase 
separation, to an elongated blob phase, to a striped phase, and finally to checkerboard-like 
phase. At suitable intermediate values of C, certain patterns are strikingly labyrinthine, as 
seen, for example, in Fig. |21for p = 0.7 and C = 3, and in Figs. El and IH for p = 0.5 and 
C = 3. 

The fact that, in the stripe regimes, the stripe widths decrease with increasing C can be 
crudely understood from the competition between long-range Coulomb repulsion and short- 
range attraction. The latter is characterized by an energy scale J, while the former has an 
energy scale C/r for spins separated by a distance r. The two are of comparable strength 
for separations r C / J . This equality occurs for larger r as C decreases - hence, to balance 
the two energies, one might expect an increasing width as C becomes smaller. 

In Table m we show the diagonal components ""y^x and of the helicity modulus tensor 7 
for the low-T spin configurations of Fig. Q, as obtained using Eq. ^ introduced above and 
discussed further below. As may be seen, the checkerboard pattern due to strong Coulomb 
repulsion dramatically suppresses 7(T) aX p = 0.5. Since 7(T) is proportional to the super- 
fluid density Us these results show that a checkerboard pattern, as expected intuitively, 
strongly suppress the superfiuid density. On the other hand, a similarly strong Coulomb 
repulsion at higher spin concentrations (p = 0.7 and 0.8) suppresses 7(T) only weakly, if 
at all. This behavior loosely resembles what is seen in doped antiferromagnets, such as the 
cuprate superconductors, as a function of doping. In this analogy, we can associate the low 
spin concentration in our model with the underdoped regime of the doped antiferromagnet, 
where it is believed that strong Coulomb repulsion leads to an insulating behavior. Like- 
wise, we can associate with the high spin concentration of our model the overdoped regime 
of the doped antiferromagnets, where the Coulomb repulsion becomes less important and in 
which they show more metallic behavior 3^. (Of course, there is no antiferromagnetism in 
our model, though there is charge ordering.) 

All of the patterns of Figs. 001 (stripes, checkerboards, and labyrinths) have been ob- 
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c 


p = 0.5 


p = 0.7 


p = 0.8 


7 


(0.00, 0.00 ) 


(0.30,0.30 ) 


(0.52,0.51 ) 


3 


(0.00, 0.07) 


(0.37,0.39 ) 


(0.49,0.54 ) 


1 


(0.32, 0.00 ) 


(0.18,0.00 ) 


(0.42,0.59 ) 


0.1 


(0.26, 0.00 ) 


(0.29,0.32 ) 


(0.37,0.42 ) 



TABLE I: Helicity moduli {"fxx^lyy) for a system described by Hamiltonian using Tc = 3a, at 
T = 0.025. C is the strength of the repulsion, and p is the spin concentration. Some snapshots of 
the system are shown in Fig. |3J 



served, both in experiments and in simulations, in systems with competing long and short 
range interactions. For example, labyrinthine structures have been observed in experiments 
on magnetic garnets 40|, and they have been obtained in simulations of spin-1 Ising Hamilto- 



„,a„s w,th co,„pet,„g long-range and short-range .nteraefons ^ Striped magnet.e pha.es 

have been observed experimentally in ferrimagnetic garnet films p|, and obtained in simu- 
lations of spin-1/2 Ising models with a long-ranged dipolar interaction !15|. Checkerboard 
patterns in the low-T electronic structure of Ca2_^Na^Cu02Cl2 have been experimentally 
observed j2^, and also obtained numerically in simulations of a classical spin-1 lattice gas 
model with short-range ferromagnetic coupling and long-range antiferromagnetic Coulomb 
interactions [IS]. There have also been numerical studies jl^ of the stripe melting transition 
in systems governed by Ising-1/2 models with short ranged ferromagnetic and long ranged 
antiferromagnetic couplings. The present model differs from all of these in having an XY, 
rather than an Ising, attractive interaction between the spins. In order to unambiguously 
distinguish our model from all these others, it would be desirable to carry out simulations 
on our model I but with an un-truncated Coulomb repulsion. 



B. Model II 



1. Numerical Results 



In our studies of Model II, we have emphasized the T-dependence of various quantities, 
in particular 7(7), as well as Si{q,T), and S2{q,T) for special values of q. We will also 
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present snapshots of spin configurations at various T's. 

Fig. shows of snapshots of the spin configurations for p = 0.5 at different T for Model 
II with A = 0, B = 0.1 (weak second-nearest-neighbor repulsion). At T = 1.2, the spins are 
already tending to form clusters of various sizes and shapes, a tendency which is clearer at 
T = 0.5. At this latter T, even though the spins form a connected path in the horizontal 
(x) direction, 'yxx{T) (shown in Fig|33) remains small, indicating no phase coherence. By 
T = 0.3, the system has phase separated into two large domains, made up of spins and 
vacancies, respectively. This transition to a phase separated configuration is signaled by 
non-zero values of Si{q,T) (Fig. Efc) for T < 0.3 and by a peak in S2{q,T) (Fig. Eli) near 
T = 0.3, with q = 2Tiy/Nya. The transition to the phase- coherent state, signaled by the 
finite value of 7(T), and that to the phase separated state, seem to occur at the same T. 

These results show that the annealed system has strikingly different behavior from a 
system with quenched disorder. In the latter case, p = 0.5 is below the critical percolation 
threshold pc — 0.59 for 2D site diluted square lattice 4l|. Hence, there is no phase coherence 
at any T in the quenched case. However, for annealed disorder, the phase separation, which 
occurs when the Coulomb repulsion is weak enough, leads to an infinite percolating cluster 
of spins, and hence to phase coherence in one of the two principal directions, even at a spin 
concentration below pc- Although similar phase separation occurs in Model I, it does not 
appear to lead to an infinite cluster of spins for p < 0.5, and hence, there is no T = phase 
coherence in that case. 

Fig. ini shows p = 0.5 for the same model but with a stronger second-neighbor repulsion 
{A = 0] B = 1.0). Fig. ink shows that at T = 1.2 and especially at T = 0.6, there is a 
tendency for stripe formation, though the stripes remain of finite length. For T = 0.55 
the system transforms to a mostly unidirectional striped phase. This transition is clearly 
signaled by a finite value in Si{q,T) for T < 0.55 (Fig. IHfc), and by a peak in S2{q,T) 
(Fig. inii) near T = 0.55, for q = nx/a. At T = 0.3, the system forms an ordered array of 
stripes in the vertical {y) direction, manifested in the fact that 5*1 (g, T) = 1 at and below 
this T. In contrast to the case of phase separation obtained for weak repulsion in Fig. (jSl), for 
which the phase coherence has a sharp onset temperature which coincides with the transition 
to a phase separated state, the phase-coherent state in the present case has a more gradual 
onset with decreasing T, and occurs at a lower T than the transition to a striped phase. 
Also, the coherence transition is observable in only one of the two principal directions, as 
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expected since phase coherence is geometrically impossible in the direction perpendicular to 
the stripes. 

We turn now to results for Model II with nearest-neighbor repulsion only. Fig. [7| shows 
results for A = 1.5, B = 0.0, and p = 0.5. As in the cases described above, the system goes 
from a disordered phase at high T, to an ordered low-T phase. However, Fig. [Tjp shows that 
the system is prevented at all T from undergoing a transition with a non-zero 7. This can be 
understood on simple geometrical grounds. Fig. [7^ shows that the strong nearest-neighbor 
repulsion decreases the number of nearest-neighbor spins at low T. Since such nearest- 
neighbors are essential in for an attractive interactions between the nearest-neighbor XY 
spins, this reduction insures that 7(T) ~ 0, as we observe numerically. The transition to a 
checkerboard-like spin pattern is indicated by the finite value of Si{q,T) for T < 0.65 and 
the peak in S2{q,T) near T = 0.65, for q = (vr/a, vr/a). 

We now present results for Model II with p = 0.8. Since p exceeds the site percolation 
threshold of pc — 0.59, 7(7) should become nonzero at sufficiently low T. Fig. |H1 shows 
results for weak nearest-neighbor repulsion {A = 0.5; B = 0.0), and p = 0.8. Fig. IHfb) 
shows that '~fxx{T) ~ '~fyy{T) for all T. We denote the helicity modulus in this isotropic 
regime simply as 7(T). 7(T) increases smoothly with decreasing T to T ~ 0.25, and, in this 
regime, coincides with 7(T) for quenched disorder at same spin concentration (continuous 
line with no symbols in Fig. (HJo). However, the quenched and annealed results differ for 
T < 0.25. 7(T) for the annealed case changes slope near T = 0.25 because the spins phase- 
separate. This phase separation is clearly visible in Fig. |HK, which shows snapshots of the 
spin configurations for several values of T. It is also signaled by the finite value of Si{q,T) 
for T < 0.2, and a peak in S2{q,T) around T = 0.2 (Figs. |Ht and|Hli), for q = 271x1 {N,j.a), 
and q = 2ny/{Nya). 

Fig. IHl shows the analogous results for stronger nearest-neighbor repulsion (A = 0, B = 1, 
and p = 0.8). In this case, the vacancies tend to cluster into stripes, not blobs - which become 
longer as T is decreased from 1.2 to 0.4; below about T = 0.2, the stripe pattern becomes 
anisotropic. The helicity moduli 'yxx{T) and 'yyy{T) become nonzero around T = 0.7 but 
remain isotropic (i. e., nearly equal) down to around T = 0.2. For T < 0.2, the anisotropic 
stripe pattern leads to a dramatic anisotropy in the helicity moduli: 'jxx{T) falls to 0, but 
'~fyy{T) increases. This behavior has an obvious geometrical explanation in the long stripes 
parallel to y, which inhibit phase coherence in the x direction. At the lowest temperature 
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of T = 0.025, there are several horizontal stripes spanning the sample, which cause '~)xx to 
vanish, while 7^,^ attains a value well above 7(T) in the case of quenched disorder. The 
transition to a striped phase is observed in the onset of a finite value of Si{q,T), as well 
as in the peak in S2{q,T), around T = 0.2, for q = (7r/a,0). In contrast to the case of 
weak nearest-neighbor repulsion shown in Fig. |S1 where 7(T) closely resembles that of the 
system with quenched disorder for all T > 0.25, the effect of strong second nearest-neighbor 
repulsion here leads to a non-monotonic 7(T): 7(T) is largest at intermediate T, and small 
at high T, where it is destroyed by thermal fluctuations, and at low T, where it is frustrated 
by the formation of stripes. By contrast, at intermediate T, short stripes have formed but 
since they are short and randomly oriented, they are insufficient to prevent a finite 7(T). 

Finally, we show the results for strong nearest-neighbor repulsion (A = 1.5, B = 0, p = 0.8) 
in Fig. ^1 Similar As in the case of nearest-neighbor repulsion at p = 0.5 shown in Fig. Q 
where the repulsion prevented phase ordering at all T, "^xxiT) and lyyiT) are substantially 
reduced compared to the corresponding quenched values, and for the same reason: the 
nearest-neighbor repulsion tends to decrease the number of spin nearest-neighbors, which 
in turn decreases the tendency of the system towards phase coherence. The Figure also 
makes apparent that an " antiferromagnetic" (checkerboard) pattern of rij's is emerging at 
low T, which competes with the XY transition. But at this p, in contrast to p = 0.5, clumps 
of checkerboard-ordered regions of zero helicity modulus can coexist with regions of finite 
helicity modulus, leading to a nonzero global value of 7(T). 

2. Analytical Results 

At T = 0, it is possible, by a simple comparison of energies, to calculate analytically the 
critical values and Be at which the system changes from phase-separated to checkerboard 
or striped order. In the phase-separated state, all the spins are contained in clusters in which 
all sites are occupied. Thus, the energy per spin in this state (assuming i? = 0) is simply 

E^' = -2 J + 2A, (8) 

where we have used the fact that each spin has two nearest-neighbor spins. For the checker- 
board ground state (taking B = 0), the ground-state energy per spin is simply 

^check ^ Q_ f^Q^ 
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The critical value of Ac is just that value of A where the two energies are equal, i. e., 
Ac = J. For A > Ac, the ground state is checkerboard; for A < Ac, it is phase separated. 
Our simulations agree with this analytical prediction. 

To calculate Be, we assume A = and compare the energies of the two spin arrangements 
at T = 0, using a simple bond-counting argument. In the phase-separated state, at T = 0, 
the ground state consists of large blobs of spins and vacancies. Disregarding the surface 
energies, we find that the energy per spin is 



since each spin has two nearest and two second-nearest spin neighbors. In the ground state 
of the striped phase (at any concentration), we assume that all the spins are contained in 
clusters consisting of alternating stripes of spins and vacancies. Since each spin in such 
a cluster has two nearest-neighbor spins and no second-nearest-neighbor spins, the ground 
state energy per spin in this phase is simply 



Be is obtained by setting these two energies equal, which gives Be = 0.5 J. For B < Be = 
0.5 J, the ground state is phase-separated, whereas for B > Be, the ground state is striped, 
independent of p. 

We have verified this prediction numerically by varying B at fixed T and several values 
of p. In agreement with eqs. (fTn|) and (jllj) . we find that the system, for any p, in the limit 
of low T, always phase separates if B < Be and forms stripes if B > Be- 

C. Low- Temperature Helicity Modulus 

For all the low-T configurations shown in Figs. (j3))- ()10p . we have compared the 7(T)'s 
obtained from Monte Carlo simulations at low T to those obtained from eq. (|7j). The config- 
urations are extracted from the snapshots at T = 0.025, and the effective conductances are 
calculated by numerically solving the system of linear equations obtained by application of 
Kirchhoff 's equations to each of the nodes in the network. To minimize finite size effects, we 
used periodic boundary conditions in the direction perpendicular to that for which we cal- 
culated the conductances. The diagonalization of the resulting matrix was carried out using 
Mathematica's built-in function 'Solve', which uses the Gaussian elimination method j^]. 



EP' = -2 J + 2B, 



(10) 




(11) 
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In Table ITTl we show the conductances for these networks, in both the x and y directions, 
as well as the corresponding values of the helicity moduli at T = 0.025. Evidently, Eq. ((7j) is 
well satisfied for the parameters considered. Where there are discrepancies, we believe that 
the source of the error is primarily the Monte Carlo simulations, since at very low T, most 
of the attempted Monte Carlo moves are rejected and the phase space sampling may be 
insufficient to give an accurate equilibrium average. We conclude that our MC simulations 
are, in general, quite accurately converged, and also that the mapping proposed in Ref. 
is well obeyed for this rather extensive series of models. 



p 


A 


B 


7x:c(T=0.025,p) 


ge,xx{p) 


7j,y(T=0.025,p) 


ae,yy{p) 


■yxx(T=0.025,p=l) 


ge,xx{p='i-) 


■yyy(T=0.025,p=l) 


9e,yy(P=^) 


0.5 


0.0 


0.1 


0.50 


0.48 


-0.01 


0.00 


0.5 


0.0 


1.0 


0.00 


0.00 


0.49 


0.48 


0.5 


0.5 


0.0 


0.50 


0.48 


0.01 


0.00 


0.5 


1.5 


0.0 


0.00 


0.00 


0.00 


0.00 


0.8 


0.0 


0.1 


0.62 


0.58 


0.67 


0.66 


0.8 


0.0 


1.0 


0.03 


0.00 


0.77 


0.76 


0.8 


0.5 


0.0 


0.64 


0.65 


0.64 


0.62 


0.8 


1.5 


0.0 


0.26 


0.17 


0.36 


0.48 



TABLE II: Comparison of the components ^xx{T = 0.025, p) and JyyiT = 0.025, p) for a site- 
diluted XY-model, as obtained by Monte Carlo simulation of Model II, to the effective conductances 
ge,xx and ge,yy of an associated conductance network, constructed as described in the text, p is 
the spin concentration while A (B) is the strength of the repulsion between first (second) nearest- 
neighbor spins. 



V. DISCUSSION 

We have studied a diluted XY model with annealed disorder and an additional spin-spin 
repulsion. We considered two types of repulsion: (i) a screened Coulomb interaction between 
spins, with a finite-separation cutoff, and (ii) a short-range repulsion. In the first case, we 
have calculated the types of minimum-energy configurations found at low T, for a variety of 
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parameter choices. For the second model, we have considered the system at finite and very 
low T. 

For the case of Coulomb repulsion (Model I), we find that, as the repulsion strength 
C increased, the system traverses a series of ordered phases in a characteristic sequence: 
first large blobs (corresponding to a phase-separated state), then horizontal, vertical and 
diagonal but straight stripes, then tortuous stripes, and finally checkerboard-like patterns. 
These patterns are strikingly independent of the model details, such how the Coulomb 
interaction is truncated. 

For Model II, the low-T spin configuration once again depends on the relative strength of 
the attractive XY interaction and the nearest-neighbor or second-nearest-neighbor repulsion. 
The ground state is always phase-separated for weak enough repulsion, but becomes either 
checkerboard or stripe-like for stronger nearest-neighbor or second-nearest-neighbor repul- 
sion. For some concentrations, we see evidence of a phase separation between a checkerboard 
or striped ordered region, and a region with no vacancies. As T is increased, the spatial 
ordering of spins, whether phase-separated, striped, or checkerboard, eventually disappears 
in favor of a homogeneous isotropic phase with only limited short-range order. The long but 
tortuous stripe regime sometimes seen in Model I at low T appears not to occur with Model 
II, probably because the repulsion is only short-range. 

Corresponding to these T-dependent spin distributions in Model II, we have seen char- 
acteristic behavior in the helicity modulus 7. For sufficiently weak repulsion, the system 
phase separates at low T. In this case, for p > pc, the two diagonal components of 7 are 
approximately equal, and substantially larger than for quenched disorder. For p < pc, phase 
separation leads to a non-zero helicity modulus in only one of the two principal directions. 
For systems with annealed disorder and strong second-nearest-neighbors repulsion, the for- 
mation of short stripes leads, at intermediate temperatures, to an increase of the helicity 
modulus in both the x and y directions. As T is reduced further, the stripes become longer 
but remain randomly oriented, leading to a reduction of the helicity modulus in both direc- 
tions. Finally, as T ^ 0, the stripes choose a preferred direction, and the helicity modulus 
becomes anisotropic, becoming large in the direction parallel to the stripes and very small 
perpendicular to them. For strong nearest-neighbor repulsion, the helicity modulus is always 
smaller than in systems with quenched disorder. This decrease is due to the reduction in 
number of nearest-neighbor spins by the repulsive interaction. In all our calculations, we 
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find good agreement between the low temperature helicity modulus as obtained from Monte 
Carlo and that inferred from the conductance of an associated conductance network. 

Finally, we comment on the original motivation for this work, which was to shed some 
light on the interplay between inhomogeneities, stripe and checkerboard order, and superfiuid 
density in the underdoped cuprate superconductors. Obviously, the present model is far too 
crude to represent all the subtleties of that system. In particular, it omits quantum effects 
arising from the non-commutativity of number and phase variablesj43|. But many of the 
phenomena reported in the cuprates (small superfiuid density, frustrated phase separation, 
and coexistence of superconductivity with other types of order, such as checkerboard or stripe 
formation) occur in our model. Thus a suitably refined version of the present model might 
provide insight into the interplay between superconductivity and other collective phenomena 
in the cuprate superconductors. 
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FIG. 1: Representative snapshots of the spin configurations at low temperatures (T = 0.025) 
for model I at various points in the p-C phase diagram. We show results for a system with an 
unscreened Coulomb repulsion truncated at r = 15a. C is the strength of the Coulomb repulsion 
and p is the spin concentration. A white (black) square is an occupied (vacant) site. Each snapshot 
was obtained by annealing the system from a T = max[2C, 1.2] as described in the text. 
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FIG. 2: Same as Fig. 1, but for a screened Coulomb repulsion with rc = 7a. 
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FIG. 3: Same as Fig. 1, but for a screened Coulomb repulsion with rc = 3a. 
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FIG. 5: Some results for Model II, with p = 0.5, and weak second nearest neighbor repulsion 
B = 0.1. (a) Representative snapshots of spin configurations at different temperatures T. White 
(black) squares are occupied (vacant) sites, (b) Diagonal components 'fxx and 'fyy of the helicity 
modulus tensor, exhibiting a coherence transition, (c) and (d): Order parameters sensitive to phase 
separation, with Qx = 217 /N^a and qy = 2'iT/Nya. 
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FIG. 6: Results for Model II, with p = 0.5, and strong second nearest neighbor repulsion {B = 1.0). 

(a) Representative snapshots of spin configurations at different temperatures T. White (black) 
squares are occupied (vacant) sites, (b) Diagonal components ^xx and 'jyy of the helicity modulus 
tensor, transition, (c) and (d): Order parameters sensitive to stripe formation, with = 7r/a, 
Qy = tt/q. 
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FIG. 7: Results for Model II, with p = 0.5, and a strong nearest-neighbor repulsion {A = 1.5). 

(a) Representative snapshots of spin configurations at different T. White (black) squares are 
occupied (vacant) sites, (b) Diagonal components 'jxx and 'jyy of helicity modulus tensor, (c) and 
(d) Order parameters sensitive to checkerboard order, with q = ^{x + y)/a. 
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FIG. 8: Results for Model II, with p = 0.8, and weak nearest neighbor repulsion (^4 = 0.5). 

(a) Spin snapshots at various temperatures T. White (black) squares are occupied (vacant) sites. 

(b) Diagonal components of helicity modulus tensor, (c) and (d) Order parameters sensitive to 
phase separation, with = 2TT /N^a and Qy = 2TT /Nya. 
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FIG. 9: Results for Model II, with p = 0.8, and strong second nearest neighbor repulsion {B = 1.0). 

(a) Spin snapshots at different temperatures T. White (black) squares are occupied (vacant) sites. 

(b) Diagonal components of helicity modulus tensor, (c) and (d): Order parameters sensitive to 
stripe formation, with = 7r/a, qy = n/a. 
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FIG. 10: Model II results, withp = 0.8, and a strong nearest-neighbor repulsion {A = 1.5). (a) Spin 
snapshots at different various temperatures T. White (black) squares are occupied (vacant) sites, 
(b) Diagonal components of the helicity modulus tensor, (c) and (d): Order parameters sensitive 
to checkerboard order, with q = Tr{x + y)/a. Note the apparent phase separation at low T into a 
checkerboard phase and a fully occupied spin phase. 
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